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ABSTRACT 



An existing computer program for the geometrically non- 
linear analysis of arbitrarily loaded shells of revolution 
was modified to incorporate initial geometric imperfections 
into the buckling solution of an axially loaded circular 
cylindrical shell. Several different boundary conditions 
were considered, and the predicted buckling loads were 
examined to determine the significance of the initial imper- 
fections and the boundary conditions for the shell stability 
problem. The computed buckling loads were compared to those 
found experimentally by Arbocz and Babcock at the California 
Institute of Technology. The computed results were 
generally higher than the experimental results. It was 
indicated that imperfection sensitivity is dependent on the 
boundary conditions of the cylinder. 
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I . INTRODUCTION 



The study of the instability of thin , isotropic , circular 
cylindrical shells has long been troubled with a lack of 
correlation between the theoretically predicted buckling 
loads and the experimental results; the experimental values 
of buckling are generally much lower than the theoretical 
solutions. This lack of agreement has been attributed to 
imperfection sensitivity, i.e., inherent imperfections 
in the shape of the shell which are unaccounted for in the 
classical buckling theory trigger the buckling modes pre- 
maturely, and to the use of incorrect boundary conditions 
in the anaiysrs , x.e. , che boundary conditions used xn me 
analysis are not equivalent to those of the test specimen [1]. 

The recent development of several sophisticated digital 
computer programs for the geometrically nonlinear analysis of 
shell structures has made it possible to include both the 
initial imperfections and arbitrary boundary conditions in 
the analysis. For example, a digital computer program for 
the geometrical nonlinear analysis of arbitrarily loaded 
shells of revolution has been developed by Ball [2]. A 
modification has recently been made to this program to 
include the initial imperfections in the analysis. However, 
this modification does not eliminate all of the previous 
analytical difficulties because in order to predict the 
buckling load of a particular shell, the imperfect shape of 
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that shell must be known, and the actual boundary conditions 
of the specimen must be known. With regard to the first 
problem, Arbocz and Babcock [3] have developed a means of 
mapping the topography of cylindrical shells and can thus 
describe the deviations from a perfect circular cylinder. 

The objective of this thesis is to use the imperfection data 
for one of the cylinders presented in Ref. [3] and Ball's 
modified computer program to obtain stability solutions for 
several boundary conditions. The buckling loads are compared 
with the experimental results presented in Ref. [3] and are 
compared with each other to determine the effect of the 
various boundary conditions. 
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II. DESCRIPTION OF THE COMPUTER PROGRAM 



The computer program used in this study is capable of 
analyzing thin isotropic shells under the following condi- 
tions : 



1) The geometric and material properties must be axi- 
symmetric, but may vary along a meridian. 

2) The initial imperfections, pressures, and tempera- 
tures must be symmetric about a meridional plane. 

3) The boundaries may be free, fixed, or elastically 
restrained. 

The governing partial differential equations are based 



on Sanders ' nonlinear shell theory for t 
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small strains and moderately small rotations (Appendix A) . 
Inplane and normal forces are accounted for but rotary 
inertial forces are neglected. The set of partial differ- 
ential equations is reduced to an infinite number of sets 
of four second-order differential equations in the meridional 
and time coordinates by expansion of the dependent variables 
in a sine or cosine series in the argument n9 where 9 is 
the circumferential coordinate and n is the Fourier index. 
Nonlinear coupling terms are treated as pseudo loads and 
trigonometric identities are used to uncouple the sets of 
equations. Derivatives with respect to the meridional 
coordinates are replaced by the conventional finite differ- 
ence approximation. This leads to sets of algebraic 
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equations in the four dependent variables , 

and where U, V, and W are the displacements in the 

meridional, circumferential, and normal directions respec- 
tively, and M is the meridional bending moment per unit 
length. At each load or time step, an estimate of the 
solution is obtained by extropolaticn from solutions at 
the previous steps. The sets of equations are repeatedly 
solved using Gaussian elimination, and the pseudo loads 
are recomputed until the solution converges. If convergence 
is not achieved in a specified number of iterations, the 
load is reduced by a factor of five. If, after a specified 
number of load reductions, the solution fails to converge 
the problem is terminated. If the load-displacement 



behavior of the shell is oi 



he 









n i" 






is presumed to have occurred at the final load. 

The original program described in [2] was not designed 
to solve shell stability problems where initial imperfections 
had to be accounted for. However, it soon became apparent 
that such a capability was very desirable and the program 
was modified accordingly. The equations and FORTRAN state- 
ments that include the imperfection terms are presented in 
Appendix A. 
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III. BOUNDARY CONDITIONS 



In the analysis, the boundary conditions for the cylin- 



der are on N or D, N „ or V, Q or W, and 3W/3s or M , 

S So S S 

where N g and N g g are membrane forces per unit length, and 
Q s is the transverse force per unit length. The boundary 
conditions are represented by the matrix equation 

r 



[fi] 



N 



N 



S0 



"\ 

U 



Q s 

3W 

3s 



>+ [A] 



V 

w 

M 

s 



= U) 



where [ft] and [A] are 4x4 matrices and {£} is a 1 x 4 
column matrix. There is one such equation for each boundary 
of the shell. Note that the analysis assumes that the 
boundary conditions are the same for all values of n , where 
n is the mode number in the circumferential direction. 

The testing configuration used in Ref. [3] for the 
axially loaded shell is shown in Fig. 1. The ends of the 
shell are attached to an end ring and a load cell with 
Cerrelow, a low melting point alloy. The ends of the shell 
are clamped such that 



W = 3W/3s = v = 0 



( 1 ) 



The axial load on the test specimen was applied by means 
of four screws located at the corners of the square end 
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FIG. I CYLIMJRICAL SHELL TESTING CONFIGURATION 
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plates. At each load increment, the screws were adjusted 
to give as uniform a circumferential load distribution as 
possible . 

In order to investigate the problem completely, several 
boundary conditions were formulated for the computer program. 
These formulations follow below. 

If the axial load is axisymmetric , then on each boundary 



N 



(0) = 

s 



N 



o 




n 1,2,... 



and 



(2) 



U (n) fl 0 



n = 0,1,2,... 



This condition will be referred to as boundary condition A. 
Note in Eq. (2) that the condition on N g is a function of n. 
This necessitates changes in the program to account for 
different values of the first element of 1 for n = 0 and 
n j- 0. These changes are given in Appendix B. 

Another possibility is 

% . ( 0 ) .. (n) « no 

SOS 

at one end of the cylinder, and 

U (n) = 0 , 0 n = 0,1,2,... 

at the other end. This condition will be referred to as 
boundary condition B. 

A third interpretation of the test set-up consists of 
rigid end plates incrementally shortening the cylinder. 

This is represented at the boundaries by 
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n 



1,2 



, • • • 



o 



at one end, and 




n 



0 , 1,2 



, • « • 



at the other end. This condition will be referred to as 
boundary condition C. 

For purposes of comparison, a simply supported cylinder 
will also be analyzed using boundary conditions A, B, and 
C. For this case, equation (1) becomes 



W = M = V = 0 



s 



(3) 
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IV. INITIAL IMPERFECTIONS 



Arbocz and Babcock's topography mapping procedure is 
described in detail in Ref. [3] . The final imperfection 
data consisted of deflections from a perfectly circular 
cylinder, denoted here by W, around fifteen circumferences 
spaced 0.5 in. apart. The first and the last circumferen- 
tial mappings occurred at a distance of 0.125 in. from the 
ends of the shell. Forty-nine data points were taken around 
each circumference with the first and the last being the 
same point on one arbitrarily chosen meridian. 

As noted previously, Ball's program requires a meridional 
plane aoour which the imperfect cylinder is symmetric. This 
requirement is not generally met by imperfect cylinders; 
thus symmetry must be created in some optimum way . An 
investigation to determine the degree of symmetry was con- 
ducted by generating a cosine Fourier series for each 
circumference using every meridian of data as a base 
meridian. In all cases, the phase angles for each circum- 
ference varied considerably, thus eliminating the possibility 
that the imperfections were symmetric. Consequently, a 
least squares method was used to determine the best meridian 
of data points about which to create symmetry. This best 
meridian was used as the origin to form the series 

W = E W ' n cos (n9 ) 
n=0 

using the data points from 0 = 0 to 0 = it. 
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An attempt was made to create a finer meridional mesh 
of the imperfection data by passing a general polynomial 
through the fifteen known coefficients for each value of n 
in order to interpolate imperfection coefficients at inter- 
mediate stations. However, a polynomial of sufficient 
degree to pass through all of the stations resulted in 
unreliable coefficients at intermediate stations. Con- 
versely, a least squares fit using a lower order polynomial 
resulted in a reasonable coefficient representation, but 
the coefficients at the original fifteen stations were 
compromised to an extent that depended on the degree of 
the polynomial. Similar problems were encountered using 
a Fourier series representation. As a result of the dif- 
ficulty in establishing coefficients at intermediate stations 
the fifteen original locations were selected as the stations 
for the numerical shell analysis and the conventional central 
finite difference scheme w as used to obtain the slope of 
the imperfection shape away from the boundaries. At the 
boundaries the simplest forward and backward differencing 
schemes were used. 
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V. AN EIGENVALUE SOLUTION 



To provide a check on the validity of the solution 
from the computer program, buckling of a simply supported 
perfectly circular cylindrical shell was formulated as an 
eigenvalue problem and solved analytically. The set of 
uncoupled field equations can be given in the form 

[E (n) ] (Z (n) "}+ (F (n) ] {Z (n) '} + (G (n) ] {Z (n) > = <e (n) } (4) 



where 



Z = {U (n ^, V (n) , W (n) , M^} 



(n ) 

The equations for the values of the elements in the [E ' ] , 

[F^], [G^], and {e^ n ^} matrices are given in Appendix C. 

For n = 0, the linear membrane solution N^ = -N is 
' s o 

assumed. Thus the buckling modes of the shell for the 
boundary conditions given in equations (2) (boundary condi- 
tion A) and (3) are 



V (n) 

W (n) 

M (n) 

s 



cos nurs/L 
C 2 sin mirs/L 
sin miTs/L 
sin imrs/L 



(5) 



where n ^ 0, are arbitrary constants, m denotes the 

axial mode number and L is the length of the cylinder. 
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Substituting equations (5) into equations (4) results in an 
eigenvalue problem of the form 

[A] {C} = -N q [B]{C} 

where [A] and [B] are 4x4 matrices defined in Appendix C 
and 



{C} = {C 1 ,C 2 rC 3 , C 4 } T . 

Because B is singular, the inverse formulation 

p- [A] {C} = - [B] {c} (6) 

o 



is used. 
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VI. SHELL DESCRIPTION AND INPUT DATA 



The shell used in this study was cylinder A-8 of Ref. 3 
Its dimensions and properties are: 



Radius, r,R g = 4.0 in. 

8.25 in . 



Nominal 

length 



Wall 

thickness ,h 

Young 1 s 
modulus , E 

Poisson' s 
ratio , v 



= 0.00464 in. 



= 15.2 x 10° lb/in. 



= 0.333 



Si v harmonics of the initial imperfection data were used 
as input to the program. They are n = 0, 9, 10, 11, 12, and 
13. The selection of these particular harmonics was based 
on the eigenvalue solution for the simply supported shell 
that predicated a minimum buckling load when m = 1 and n = 9, 
and the experimental observation in Ref. [3] that n = 13 
appeared to be the critical mode where m and n are the mode 
shapes in the meridional and circumferential directions 
respectively. The number of meridional stations used was 
fifteen because no reliable interpolation of the initial 
imperfections could be developed, as discussed above. The 
length of the shell, L,was taken as 8 in. with the first and 
the last data point stations assumed to be the stations at 
the edges of the shell. Convergence criterion was taken 
as 0.01. 
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As presented in Appendix A, the imperfection input data 



consists of the rotations <f> and <f>^ n ^ where 




/awM 

J* 



and 



(n 

e 



ao 

o 

rE o 




in which o is a reference stress level, and E is a 
o o 

reference elastic modulus, a is a reference length, and k 
denotes the station number. The data are introduced into 
the program through the subroutine IMPERF. The FORTRAN 
statements for IMPERF are given in Appendix D. The imper- 
fection data was provided to the author by the writers of 
Ref. [3] in the form of a deck of cards with the data W/h. 

In order to obtain a buckling load that was a very 
close approximation to the bifurcation lead or tne eigen- 
value problem, the imperfections were reduced by a factor 
10 This technique has been used successfully with this 

program by Stillwell [4] who used very small asymmetries 
in the load as the triggering mechanism. 
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VII. RESULTS AND DISCUSSION 



The solution to the eigenvalue problem given by equa- 
tion (6) resulted in a minimum buckling load for the simply 
supported cylinder of 53.02 lb/in. when m = 1 and n = 9. 

The solution from the computer program for the minute imper- 
fections resulted in a buckling load of 54.0 lb/in for 
boundary condition A, simply supported edges, and displace- 
ment in harmonics n = 0 and n = 9. Some difference between 
the two results is to be expected because the analytical 
solution neglects the pre-buckling deformation due to edge 
constraint, and the numerical solution used only fifteen 
s tations . 

The computer results for the simply supported cylinder 
with boundary conditions A, B, and C, with both minute and 
full imperfections,* are presented in Table I. These results 
show a relationship between the imperfection sensitivity of 
the shell and the boundary conditions. For example, boundary 
condition A shows a 15% reduction in the load carrying capa- 
city when the full imperfection is introduced. Boundary 
condition B shows a 10% reduction, while condition C reduces 
the capacity by 42.7%. 

The computer results for the clamped shell are given in 
Table II. The experimental buckling load of shell A-8 was 

* -6 
minute = full imperfections x 10 
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TABLE I 





SIMPLY SUPPORTED 


CYLINDER 


Loading 

Condition 


Degree of 
Imperfection 


Buckling Load 
(lb/in) 


A 


minute 


54.0 


A 


full 


45.9 


B 


minute 


69.0 


B 


full 


62.7 


C 


minute 


78.6 


C 


full 


45.0 



TABLE II 

CLAMPED CYLINDER 


Loading 


Degree of 


Buckling Load 


Condition 


Imperfection 


(lb/ in) 


A 


minute 


7.5 


A 


full 


no 






convergence 


B 


minute 


43.87 


B 


full 


46.9 


C 


minute 


85.7 


C 


full 


49.1 
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found to be 32.87 lb/in [3]. From Table II: 

1) The buckling load for condition A with minute imper- 
fections is 7.5 lb/in.* 

2) No converged solution could be found for condition 
A with full imperfections. 

3) The buckling load for condition B is slightly higher 
with full imperfections than for the minute case. 

4) The buckling load for loading condition C is 
decreased by 42.8% when the full imperfections are 
introduced. 

Thus, the stability of the clamped shell under condi- 
tion A is drastically reduced below that of condition B and 
C. This is very puzzling and much effort was expended in 
search of an explanation. Unfortunately, none war, rouna. 

The computed result for condition C, although generally 
higher than the empirical result, shows the expected load 
reduction due to the initial imperfections [3], Some of 
the discrepancy in magnitude could be due to the use of an 
effective length that was 0.25 in. less than the actual 
length of the shell since the greater length would lower 
the computed buckling load. The remainder may be due to 
the fact that the actual boundary condition is some combin- 
ation of conditions A and C. Another consideration is that 
the symmetric representation of the actual imperfections 
may also degrade the validity of the solution. 

* 

For this boundary condition, the imperfections were 
reduced by a factor of 10“9. 
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The circumferential harmonic n = 10 was found to exhibit 
the most rapid growth. Figure 2 displays vs for 

boundary condition C. The sudden decrease in slope occur- 
ring just prior to reaching illustrates the buckling 

(non-convergence) phenomenon. 
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N (0) VS W (l0) STATION 9 
s 




FIG.2 
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VIII. CONCLUSIONS 



1. The imperfection sensitivity of a cylinder is dependent 
on the boundary conditions. 

2. The computer program yields more realistic results for 
the simply supported boundary condition than for the 
clamped condition. 

3. The initial imperfections result in a significant 
reduction in load carrying capacity in the simply 
supported case but less reduction than expected [3]. 

Further investigation is necessary to determine: 

1. Why are Lhe computed l -salts generally higher than the 
experimental results? 

2. Why is the clamped boundary condition A so sensitive to 
imperfections in the computer program? 

3. Why does loading condition B, clamped, appear to be 
strengthened by the introduction of full imperfections 
in the computer solution? 
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APPENDIX A 



This section presents the appropriate equations for 
including initial imperfections in the computer program. 

The strain-displacement relationships used in this analysis 
are 

£ = U' + W/R + ($ 2 + $ 2 )/2 

3 3 3 

£q = V'/r + r'U/R + W/R + ($ 2 + $ 2 )/2 (Al) 

e sQ = (V' + U’/r - r'V/r + $ c ,$ Q )/2 

where 

$ = -W ' + U/P. 

S 3 

=-W*/r + V/R q (A2 ) 

$ = (V' + r'V/r - U'/r)/2 

and where e , e Q , e Q represent Fourier coefficients of the 

S U So 

reference surface strains, $ , $ 0 , and $ s0 represent refer- 
ence surface rotations, and R 0 and R g are the principal 
radii of curvature. The superscript primes and dots denote 
partial differentiation with respect to s and 0 respectively. 

Consider an imperfect cylinder. Before loading, the 
imperfections appear as a displacement from the perfect 
cylinder. Note that this is a displacement with no accom- 
panying strain. Call this displacement W, and the total 
displacement from the perfect shape W* . Then 
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w* = w + w 



where W is that part of the total displacement caused by 

the applied load and is associated with strain. Replacing 

W with W* in equations (Al) and (A2) and requiring that 

e = e n = e a = 0 when U = V = W = 0, leads to 
s 9 s0 

e g = U'+W/R s + [(-W'+U/R s ) 2 - 2W' (-W+U/R ) +$ 2 ]/2 

e Q = V/r+r'U/r + W/R 0 + [(-W'/r+V/R ) 2 - 2W*/r (-W ' /r+V/R^) 
+ $ 2 ]/2 



e = [V'+U'/r -r'V/r + (-W'+U/R ) (-W'/r+V/R ) - W' 

S0 S u 

(-W* /r + V/R Q ) - W*/r(“W' + U/R ) ]/2 



Defining the imperfection rotcti 

$ = — W ' and 

s 9 



a r> H ^ c 

“ *«? ~ '9 

-W'/r 



equations (A2) become 

£ = U'+W/R + ($ 2 + 2? $ + $ 2 )/2 

s s s s s 

e 0 = V‘/r + r'U/r + W/R Q + ($ 2 + 2? 0 $ 0 + $ 2 )/2 
£ = (V ' +U* /r - r ' V/r + $ $ Q + $ + $ $„)/2 

SQ ' S0S0S0 



From Ref. [2] 



$ 



2 

s 





cos 




a 00 , s 

=£ Z 3 n cos n0 
E o n=0 s 



Redefine 3 such that 
s 
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0 00 

= 2 - E B tn) cos n 6 = •S ' 2 + 2 ? * 
E o n =0 s s s s 



°o II " , <n) 

E* S E +s cos n0 
o \n =0 



+ 2 



00 \ / 00 
Z (j)^ n ^cos n 0 | E <j>^cos n 0 
n =0 s \ n=0 s 



Similarly 



o ~ 9 

=T- Z 6 COS n 0 = $n+ 2 
E o n =0 6 0 e e 






„ . n . 

Z <PqCOS n 0 

n=l 



^ + 2 ( Z sin n ®|( ^ s ^ n n ®l 



and 



a 00 

=r^- 2 6 fl sin n 0 = $ + $ <J> + $ 

E , s 0 s 0 s 0 s 0 

o n=l 



/ ' ft 

f a 



. ' t \ 

o _ \ r> ! / °" , , « • «" x 

Z 4>^ n; cos noli Z ^ h, sin n0 

n =0 s n=l 6 







Z cfi^cos n 0 )( Z <()p n ^sin n 0 1 + 
n =0 s \n=l W 



Z -r(n) • o 1 

Z <j> ' sm n 0 , 

n=l b 



Z (j> cos n 0 | 
n =0 S 



where 



-r(n) = _ dw 



“(n) 



dK 



and 



— (n) nw^ n ^ 

*e “ -T~ 



:(n) 



, „ / / , — (n) = E W v /(a o ). 

and c, = s/a, p = r/a, and w o o 



The changes to the program are in SUBROUTINE PHIBET 
given below. 
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SUBROUTINE PHI BET ( K ) 

COMMON 

1 / I B L 1/MN M A X 
2 / 1 B L 2 / N ( 10 ) , MN IN I T 
3/IB LA/KMA X , KL 

3/IBL7/MNMAXQ,MAXD (10), MAXS( 1C) ,MAXSY< 10),IS(10,10),JS(1G,10),ID(10 
4,10 , JD( 1C , 10) , IJS( 10 ) 

6/1 BL 12 /KM A X 1 , KMA X2 , NCONV 
2/IBL13/ITRMAX, LSMAX 
8/BL6/Z<4,220) , SOE , OSE , ALOAD 
0/BL8/RI 2 OC ) ,OAM( 20") , OMT( 2uO) 

9/BL 10/PHI V ( 10 ), PH IT (101 , PHI < 10) 

3/BL1 1/OMXI (200) ,PHEE,T0,T2 
1/BL12/T0L I , TOEL 

2/BL15/NU,Ul (10) tVl ( 10) ,W1 ( 10) , V2( 1C) ,U2( 10) ,W2( 10) ,U3( 10) , V3( 10) , 
3W3( 10) 

4/BL27/BX3( 10) . BT3 ( 10) , BXT3( 10 ) ,BE3( 10) 

COMMON/S LI MP/PHIXB(200 ) , PHITB(200 ) 

COMMON /BL P HS/ PHX ( 10) ,PHT( 1C) 

CX=CMX I ( K ) 

OT=CMT ( K ) 

RRA=1»/R(K) 

GA=GAM(K) 

KP2=K+2 

DO 1 M=1 , M r.'MAXG 
EN=N( M) 

IK=KP2+( M- 1 )-KMAX2 
U3(M)=Z( 1 , IK) 

V3(M)=Z(2 , IK) 

W3(M)=Z (3, IK) 

PH I X ( M ) =-T OL I * ( V/3 ( M ) -w 1(M) )+0X*U2(M) 

PHI T ( M ) = EN*W2 ( M) t? RR A+V2 ( M ) *0T 

1 PHI (M)=(TDLI-MV3( M)-Vl(M) )+GA*V2(M) tEN*U2( M)*RRA)*.5 
I F ( ITRMAX . EG . 1 ) RETURN 

DO 60 F = 1 , MMMAXO 
KP=K+ ( M~ 1 ) *K W AX 
PHX(M)=PHIX6(KP) 

60 PHT(M)=PHI TB(KP) 

Qn Q ,v=l . MNMAX 

SMO*G. 

SMT=0. 

SMR=C * 

SMF =0 e 

IF(N(M) . EO.C ) C-0 TO 20 
MAXL=MAXS ( M ) 

IF ( MAXL* EO. C ) GO TO 2 
DO ? L= 1 , MAXL 
I=IS(L,M) 

J=JS(L,M) 

SMO=SMO+PHIX( I )*PHIX( J ) 

1 + PH X ( I ) 5 PH I X ( J ) +PHX( J )*PHIX( I ) 

SMT=SMT-PH IT ( I ) sV PHl T ( J ) 

1 - PHT(I )*PHIT (J ) -PHT( J)*PHIT( I ) 

SMR=SMR+PH I X ( I ) *PH I T ( J ) + PH I X ( J ) *PH I T ( I ) 

1 +PHX ( I )*PHIT( J ) + PHX( J)*PH!T( I )+PHIX( I ) * PHT ( J ) +PH I X ( J ) *P HT ( I ) 

3 SMF=SMF-PH I ( I ) *PH 1 ( J ) 

2 MAX L=MAXD ( M ) 

IF(MAXL.EG.C) GO TO 4 
DO 5 L = 1 » MAXL 

I = I D ( L , M ) 

J=J P( L , M ) 

SM0 = Si10+PHIX( I )*PHIX( J ) 

1 + PHX ( I ) *PH I X ( J ) + PHX ( J ) ’rPHTX ( I ) 

SMT=SMT+PHTT ( I )-PHIT( J) 

1 + PHT( I)*FHIT( J) + PHT ( J )*PHIT ( I ) 

SMR=SMR-PHIX( I )*PH I T ( J ) +PHI X ( J ) --PH I T( I) 

1 -PHX ( I ) vPH.IT ( J )+PHX(J)” ; PHIT( I ) -PH I X ( I ) *PHT ( J ) + PHI X ( J )*PHT( I ) 

5 SMF = SMF+ PM I ( I ) =>■ PH I ( J ) 

4 I F ( MAXSY ( M ) . EQ. C ) GO TO 10 
1=1 JS(M) 

SMO=SMG+PHI X( I )**2/2. 



30 



1 +PHX( I )*PHIX ( I ) 

SMT = SMT-PH I T ( I )**2/2. 

1 -PHT(I )*PH1T(I ) 

SMR=( SMR + PH1X( I )*PHIT ( I ) ) 

1 + PHX ( I ) * : PH IT ( I ) +PH IX ( I ) =4'-PHT ( I ) 

SMF=SMF-PHI ( I )**2/2. 

GO TO 10 

0 CO 21 L= 1 , f. NMAXO 
SKC=S"C + PHIX(L )**2 

1 +2.vpHX(L)--'PHIX(L) 

SMT = SMT + PHIT(L )**2 
1 +2.*PHT(L)*PHIT(L) 

1 SMF=SMF+PH I (l_ ) "=> : 2 

I F ( GT • MNNAXO ) GO TO 11 
SMG=SMC+PH1XI,’<)**2 
1 + 2 • P H X ( f * } * • P FI I X ( M ) 

1 BX3 ( P ) = Sf\0*« 5 
BT3(M)=SKT*.5 
BE3(M)*SMF*.5 
BXT3(tf)=C. 

GO TO 9 
G BX3 ( M ) =S MO 

B T 3 ( ) = S M T 

BXT3(M)=S^R*.5 
BE3 ( N ) =SMF 
9 CONTINUE 
RETURN 
END 
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APPENDIX B 



The load applied in the computer program is uniform and 
taken in the first mode only, i.e., n = 1. To accomplish 
this, the following changes were made to the program logic: 



SUBROUTINE ZANDZ 531 



DO 15 J = 1, 4 659 

* IF(M.NE.l) ELLS (J) = 0 



RETURN 

END 



SUBROUTINE FORCE (K) 960 

DO 21 J = 1, 4 081 

* IF(M.NE.l) ELLS (J) = 0 



RETURN 

END 



* = Required inserts 
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APPENDIX C 



The elements of the [A] and [B] matrices used in the 
eigenvalue solution are defined below. 
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where Q = rms/L 



and 
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APPENDIX D 



The following subroutine was used to calculate <J> ^ 
_ (n) S k 

and <J>q ; in the computer program, 
k 



SUBROUTINE IMPERF 
DIMENSION FNT (15) , DFNT (15) 

COMMON /BLIMP/PHITB (200) ,PHIXB(200) 

AO=4 

SIGA=1000. 

E0=15200000 . 

H=8./14. 

Hl=l./H 
H2=l./(2. *H) 

DO 1 1=1,6 
AJ2=I+7 

IF(I.NE.l) GO TO 8 
AJ2=0 . 

S READ (5,10) (FNT Ml , J =i , 15 ) 

DO 2 K=2 , 14 

2 DFNT (K) = (FNT (K+l) -FNT (K-l ) ) *K2 
DFNT (1) = (FNT (2) -FNT (.1) *H1 
DFNT ( 15 ) = (FNT (15) -FNT(14) ) *H1 
KR=I-1 

DO 3 Kl=l ,15 
KZ=KR*15+K1 

PHITB (KZ)= (FNT (K1 ) *AJ2*EO/SIGA) *(0.00464) 
PHIXB (KZ)= (DFNT (Kl) *EO/SIGA) *(-0.00464) 

3 CONTINUE 
1 CONTINUE 

10 FORMAT (6E12. 5) 

RETURN 

END 
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